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Abstract 

We present the extension of Frozen Density Embedding (FDE) theory to real-time 
Time Dependent Density Functional Theory (rt-TDDFT). FDE a is DET-in-DET em¬ 
bedding method that allows to partition a larger Kohn-Sham system into a set of 
smaller, coupled Kohn-Sham systems. Additional to the computational advantage, 
EDE provides physical insight into the properties of embedded systems and the cou¬ 
pling interactions between them. The extension to rt-TDDFT is done straightforwardly 
by evolving the Kohn-Sham subsystems in time simultaneously, while updating the em¬ 
bedding potential between the systems at every time step. Two main applications are 
presented; the explicit excitation energy transfer in real time between subsystems is 
demonstrated for the case of the Na4 cluster and the effect of the embedding on optical 
spectra of coupled chromophores. In particular, the importance of including the full 
dynamic response in the embedding potential is demonstrated. 
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1 Introduction 


It is the nature of a chemist to view molecules as a collection of atoms and condensed states 
as a collection of molecules, and to interpret the properties of the larger aggregates in terms 
of their smaller building stones and the interactions between them. In the held of computa¬ 
tional chemistry, this tendency was hrst expressed in a development of a range of so called 
Atoms-In-Molecule methods which are meant to be used as a post-processing 

tool in order to reveal the subtle relations between structure and properties, hidden in the 
single quantity obtained from a supermolecular calculation. An array of AIM methods based 
on space partitioning,^®^ density partitioning and wave function partitioning have been 
developed over decades and successfully applied to a range of properties as charges, ener¬ 
gies, response properties and reactivity analysis. Embeddin^®®®!!^^^ is a more recently 
developed held which can be seen as taking AIM a step further and instead of applying its 
concepts in the post processing, merge it with the electronic structure theory. In this process, 
the goal to seek deeper chemical understanding has been supplemented by the aspiration to 
increase computational efficiency. Embedding has also opened the door to combining dif¬ 
ferent electronic structure methods into a single calculation, building a bridge between the 
worlds of wave function theory and density functional theory (BET) SIllini P^W IIsllIB] 

Subsystem BET, also referred to as Frozen Bensity Embedding (FBE) theory (see Section 
12.11 for discussion), is closely related to the Hirshfeld AIM method,® as both view the total 
density at each point of space as a sum of subsystem densities 

Ns 

Ptot(r) = ^p/(r). (1) 

I 

Where the Hirshfeld method handles each atom as a separate entity, FBE partitions the 
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subsystem into molecular fragments. Though Eq. ([T]) does not seem to imply it at first sight, 
the partitioning also includes the division of space. The densities of the subsystems can 
overlap and, while they are in principle allowed to be delocalized over all space, in prac¬ 
tice the initial guess is based on a chemically sensible representation which usually leads 
to a collection of weakly overlapping densities, centered around the nuclei of the subsys¬ 
tems.!^^ The “cutting” of covalent bonds is avoided, since the description of the interaction 
between the subsystems relies on the use of approximate nonadditive kinetic energy function¬ 
als (NAKE), which, at this point, still struggle with strongly overlapping densities.!^ This 
through-space partitioning is the key to the simplihcation of the computational problem in 
FDE, which allows to take advantage of computational techniques such as nearsightedness 
and the use of localized basis sets. In this sense FDE succeeds in combining the two goals 
of deeper chemical insight and computational efficiency.!^!^!^!^!^!^ 

Most embedding methods are developed to be used for ground state calculations, though 
recently several have been extended to the realm of time dependent methods.®!^!^!^®!!^ 
The hrst extensions have been made for the linear response formalism of time-dependent 
DFT (LR-TDDFT),!2S1I23 ^ theory developed for the calculation of excited states that has 
increasingly gained in popularity in the last two decades.!^!^!^ Real-time time-dependent 
density functional theory (rt-TDDFT), on the other hand, has only relatively recently came 
into the spot light ,!2SM36J137i[38] j^iostly due to its role as a stepping stone into the world of 
nonadiabatic dynamics. In contrast to LR-TDDFT, rt-TDDFT does not rely on the linear 
response formalism and is intrinsically able to describe the full response of the system to an 
applied perturbation.^^^l With the introduction of the Kohn-Sham DFT (KS-DFT) formalism 
and the adiabatic approximation (see Section 12.21 for discussion) it reduces to a non-linear 
integration in time of the time dependent Kohn-Sham (TDKS) equations. While this can 
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be computationally demanding, as it requires very small time steps 1-2 attoseconds (as)] 
and long simulation times |]~ 10-200 femtoseconds (fs) depending on the application], it has 
the advantage of scaling as 0{N^). While this implies that linear scaling techniques can be 
applied as in the case of ground state KS, one must be cautious with introducing numerical 
errors as they can accumulate during the long simulation process. 

The integration of rt-TDDFT in embedding methods has only started recently and, until 
this work, was limited to the exact embedding method^^ The spirit of these methods 
is very different from FDE, as they require a precalculated total density, aiming at com¬ 
bining different computational methods. Here we present a new formulation of subsystem 
rt-TDDFT and show (see Section 14.21 for discussion) that combining rt-TDDFT with a 
subsystem approach can offer a view into the process of excitation energy transfer, which 
cannot be obtained straightforwardly from standard rt-TDDFT calculation.^^ Additionally, 
subsystem rt-TDDFT readily provides information about the role of coupling between the 
subsystems in the reproduction of optical absorption spectra, a question that has been raised 
before in LR-TDDFT applications.^^ 

2 Method 

2.1 Subsystem DFT 

In this section we will review the fundamentals of subsystem DFT before extending the 
formalism to rt-TDDFT. In subsystem DFT the partitioning of the total system into sub¬ 
systems is made at the density level [Eq. ([T])]. In DFT, the density is the central quantity 
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and is used to formulate the variational principle for the ground state energy of the system 


Eq = min E[p] = min (F[p] + / next(r)p(r) ) . (2) 

p^N p^N \ J J 

where, using Levy’s constrained search, p is restricted to be A^-representable, Vexti'^) is the 
external potential due to the electron-nuclei interaction and E[p] is the universal functional. 
E[p] is dehned as 

FM = Tip] + H'IpI = T.ip] + J|p] + E„[p] (3) 

where T[p] is the kinetic energy and W[p] is the electron-electron interaction. F[p] can also 
be rewritten using the noninteracting kinetic energy Ts[p]., the Coulomb energy J[p] and 
the exchange-correlation energy Exc[p] , dehned as the difference between T[p] and Ts[p] and 
W[p\ and J[p\. In the case of KS-DFT, the density p is mapped to a set of N non-interacting 
electrons represented by a single Slater determinant consisting out of the KS orbitals 0j(r) 
and is given, for closed-shell systems, by 




(4) 


This allows to rewrite the energy expression as 


Eo = min (E[p]) = min (T,[{(t)i}] + [ Veff{r)p{r)\ . 

4>^p^N ^ \ J / 


(5) 


where the noninteracting kinetic energy is now expressed exactly using the KS orbitals 


N 


T.m] = E 


i=l 


( 6 ) 
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and all electron-electron interaction terms are stacked away in an effective potential 


e//(r) = Vext{r) + - + 


5J[p] , 5E,,[p] 


5p{y) 5p{r) 


Minimizing E[p] with respect to the KS orbitals is thns achieved by 



N 

E[p]-Y^ ep 

o 



0*(r)(^i(r)dr 



0, \/k 


which leads to the one-electron KS eqnations 


( 7 ) 


( 8 ) 


-^V^ + Veff{r) 


(r) = 


(9) 


When the total density is decomposed into a snm of snbsystem densities {^/(r)} [Eq. 
O], the total energy fnnctional E[p] can be rewritten as 

E[p] = E[{pi}] =J2[ '^ext{r)pi{r)dr+J2Yl 

I 1 j 

( 10 ) 

where the external potential term is additive, the Conlomb repnlsion term is pairwise additive 
and the kinetic energy and exchange correlation energy are nonadditive. The nonadditive 
density fnnctionals can be rewritten as 



pi{r)pj{r’] 


r — r 


drdr'+Ts 


Up- 

. I 


+Ex 


, I 


K 



^ K \p,] + hf \p]-'£k [p,] j Ip,] + If”--'-' |{p,}] 


( 11 ) 


If the division into snbsystem densities is done in snch a way that each one is noninteracting 
v-representable, then they can be viewed as separate KS snbsystems, each mapped to an 
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effective KS potential and a set of subsystem KS orbitals The additive part of the 

noninteracting kinetic energy functional in Eq. flTT|) can then be expressed exactly using the 
subsystem orbitals Ts[{0f}] = while the nonadditive part is expressed 

using a kinetic energy functional Ts[{pi}]. In the spirit of KS, all other terms in Eq. flTU]) 
are gathered together in the subsystem effective potential vfyy(r) after taking a functional 
derivative to pj 


Sp{r) Sp{r) Sp{r) Spi{r) 


( 12 ) 


Minimizing the energy of the total system is then replaced by a set of coupled minimizations 


En = minmin... min (Eljorjl). 

m {4>n 


(13) 


Each of the minimizations is equivalent to 


EliPi}] 






4'f 


= 0 


(14) 


which leads to a set of coupled subsystem KS equations 




(16) 


Derivation of the subsystem KS equations depends on the following constraints for each 
subsystem /: 


1. The density is noninteracting pure-state v-representable (fs-representable). 


1. 5pj = 0, J 7^ /, to insure that = 5(r - r') and = 6ij6{t - r'). 











Due to constraint HI Eqs. fllSp are solved iteratively through “freeze-and-thaw” cycles]^ 
where the energy is minimized with respect to subsystem I while keeping the densities of 
all other subsystems frozen. At self consistency, the energy is minimized with respect to the 
variation of all subsystem densities and therefore also the total density. 

The subsystem effective potential in Eq. flT^ can be rewritten as the KS potential of the 
total system plus the nonadditive kinetic energy potential 


T// 


,tot 


T) = Veff 


(r) + 


Sfsjp] 

5p{r) 


St[pi] 

6pi{r) 


(16) 


where we use the notation Tg to emphasize that an approximate density functional is used. 
In case of the exact kinetic energy functional Ts[p] the subsystem effective KS potential 
would reduce to the exact KS potential associated with the embedded density p/(r) plus a 
constant 

^i//(r) = + const (17) 

where we used the relationshippSl 


STgjp] 

Sp{r) 


const — u^j(r) 


(18) 


In other words, in the limit of exact NAKE, FDE reproduces the exact density and KS 
orbitals of the embedded fragment. However, in practice, approximate NAKE functionals 
are used for the nonadditive kinetic energy and the subsystem KS effective potential differs 
from the the exact KS potential b})^ 


I ( ^ ^ ^ ^Tg[pi] _ 5fg[pi] 

6p{r) Sp{r) hp/(r) Spi{r) 


(19) 
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The performance of the approximate NAKE is therefore essential for the reproduction of the 
densities of the embedded fragment and, as a result, also the total supermolecular density. 

During the freeze-and-thaw (FAT) cycles, the density of subsystem I is being repeatedly 
optimized, each time in the presence of a different frozen density constraint. At each FAT 
cycle, the self consistency corresponding to a different KS potential Ugjy(r) is found. This 
process can be considerably sped up by updating the KS potential after each SCF cycle.^ 
Since the self-consistency point is well defined, both approaches lead to the same result. 

Subsystem DFT has been shown to perform well for systems where the density can be 
partitioned between non-covalently bonded systems, i.e., systems interacting through elec¬ 
trostatics, van der Waals forces and hydrogen bonds.^®®!^ Covalently bonded fragments, 
as well as systems with partial charge transfer character, remain a challenge for subsystem 
DFT until more accurate NAKE are developed.®! When applied to systems with clearly 
separated fragments, subsystem DFT offers the advantage of providing us with a chemically 
sensible approximation of the density and properties of embedded systems. The embedding 
potential, given below, represents the difference between the KS potential of the noninter¬ 
acting fragment, (i.e., evaluated with the density of the embedded fragment but not 

including the interaction with the environment), and the effective KS potential of the em¬ 
bedded subsystem As a result, it contains the full information about the interaction 

and coupling between the subsystems. 


7 / 
^emh 


(r) = vijfiv) - uLnt(r) = XI' 

j^i 


E 

j^r 


^dr' + 


nadd 


r — r 


+ 


5t: 


nadd 

s 


Spi{r) Spi{r) 


( 20 ) 
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2.2 Subsystem real time TDDFT 


The central point of the derivation of snbsystem DFT for the ground state is the variational 
principle, i.e., minimization of the total energy of the system with respect to the density of 
each subsystem. In contrast to the ground state case, there is no such minimum principle 
for the time dependent problem and therefore no unique way to derive the time evolution 
of a system. A possible route is to use the stationary action principle which states that the 
variation of the quantum mechanical action A, dehned as: 


.41*] 





H(t) 


*(*) 


( 21 ) 


should be zero. In Eq. H is the full molecular Hamiltonian that includes the time- 

dependent applied potential and \h(f) is the electronic wavefunction. By imposing the 
boundary conditions that hT(fo) = (5T(ti) = 0 one recovers the time-dependent Schrodinger 
equation 

= mm ( 22 ) 

Since the introduction of the Runge-Gross theorem,^ there has been ongoing discussion 
whether the stationary action principle can be carried over straightforwardly to TDDFT. In 
TDDFT, the (fully interacting) quantum state \h(t) is a functional of the time dependent 
density making the action a density functional as well. 


A[p\ 


'to 


dt{^p]{t) 


rti 


'^0 




. d ^ 
'di ~ 


W 


mit)) 





drp{r,t)vext{r,t) 


(23) 

(24) 
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As was shown by van Leenwen,l^l^^ the density p(r, t) which evolves in an interacting 
system under the influence of an external potential nea;t(r,t) starting from an initial state 
^[p](^o) can be reproduced in a KS noninteracting system evolving under the action of an 
effective potential Vs{r,t) starting from an initial state <h[p](to), where <h[p] is the single 
Slater determinant representing the noninteracting system. 


AM 




AM 






drp{r,t)vs{r,t) 


(25) 


The difference ^^(r, t) — vnir, t) — Vexti'^, t) defines the time dependent exchange-correlation 
potential Vxc{i^,t). The term nea;t(r,f) contains here both the nucleus-electron interaction as 
well as any applied time-dependent held in the form of 


Vappi{r,t) = vo{r) + Vi{r, t)9{t - to) (26) 

where 9{t — to) represents a step function. 

The variation of the TDDFT action is unfortunately not as simple a matter as in the 
wave function formulation of the theory. Since the action is a density functional, one needs 
to take the variation with respect to the density. The boundary conditions are, however, not 
transferable. One can still require to be zero for to but not for ti, since a variation in 

p at a time to < t < ti can and will change the quantum state <hM at time ti. The (not 
longer that) stationary action principle can then be rewritten as 

6 A[p] - i (d^M(ti) I ^^M(ti)) = ^AM - i (*^’M(ti) I ^*^’M(ti)) = o (27) 
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The problem is circumvented in the usual spirit of KS-DFT by shifting the unknown terms 
into the exchange-correlation potential, explicitly dehned as 


^ 5A^c[p] . (T[p](ti) |(5T[p](ti)) _ . (<h[p](ti) |^<l>[p](ti)) 

5p{r,t) ^ 5piy,t) * 


In practical calculation, one usually resorts to the adiabatic approximation, which simplihes 
the time-dependent exchange-correlation potential to the ground state exchange correlation 
potential evaluated at the time dependent density. 




^Exc [p] 


p=p{r,t) 


(29) 


As a consequence, the extra terms in Eq. fl25]) vanish, the density time dependence is reduced 
to instantaneous and causality is trivially fulhlled. The familiar one-electron time dependent 
KS equations (TDKS) take the form 




(30) 


where 0j(r,t) are the noninteracting KS one-electron orbitals constituting the single Slater 
determinant ‘h[p](t) and yielding the time dependent density p(r,t) 

Once a solution to the ground state coupled subsystem KS equations [Eq. flTHl) ] has been 
obtained, we can dehne an action principle for each of the subsystems I. Each subsystem 
is represented by a single Slater determinant $[p/], which is a functional of the subsystem 
density p/(r,t). Namely, 





._a 


^[p/]W^ 



drp 7 (r,t)nf(r,t), (31) 
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where, in accordance with the adiabatic approximation, Vg{r,t) is dehned as: 


vl{r,t) = Vext{r,t) + 


dJ[p] , dExM , St[p] dt[pi] 


+ 


+ 


6p{r,t) Sp{r,t) Sp{r,t) Spi{r,t) 


(32) 


and is equal at t = to to the in Eq. flTB]) at self consistency. We note that only 

t) contains an explicit time dependence while the other terms depend on time through 
the density. In the limit case of exact NAKE, the sum of the subsystem actions [Eq. fl3T]) ] 
will reduce to the total KS action [Eq. fl2S]l ]. Taking the variation of each subsystem action 
allows us to obtain the subsystem TDKS equations 


d 




dt 


(33) 


In order to reproduce the total KS action, it is imperative to integrate the subsystem 
TDKS equations simultaneously. The subsystem action, as dehned in Eq. (1^ . represents the 
coupled action, since the subsystem potential vf(r,t) depends on the total time dependent 
density p(r,t) = system. If the applied time dependent potential to the 

system is sufficiently small, one can draw a parallel to the subsystem formulation of linear 
response TDDFT, where the coupled response function of the subsystem yf = is 

connected to the coupled response functions of all the other subsystems through a Dyson- 
type equation (in simplihed notation)^*^ 

Ns 

= x'S + Y. xWux‘j ( 34 ) 

J^I 

where Kjj = As in the case of subsystem formulation of linear response TDDFT, 
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one can also define an equivalent to the uncoupled response of the subsystem 


xy = X? + X^iKjiXl (35) 

where x? is the subsystem KS response function. This is achieved by dehning the uncoupled 
subsystem action, where the density of the other subsystems are kept frozen and Eq. fl33l) is 
only evolved in time for the subsystem in question. As we will show in Section ITTl examining 
the difference in subsystem properties between the coupled calculations and the uncoupled 
calculation reveals hrst hand information on the excitation coupling between the subsystems. 
It is also important to note that the embedding potential is never fully static also in the 
uncoupled calculations, as it depends on the total electron density which varies in time even 
when the other subsystems are kept frozen. 

3 Implementation details 

The subsystem rt-TDDFT method was implemented into the Quantum Espresso package 
as an extension of the periodic subsystem DFT®! using plane waves, the details of which can 
be found in Ref. [23] In this work we only report calculations involving molecular systems 
using the T point sampling for the Brillouin zone. All calculations have been performed 
using ultrasoft pseudopotentials.l^l^l^ Ultrasoft pseudopotentials can be seen as a special 
case of the projector augmented wave (PAW) methocP^ which reformulates the KS equations 
in terms of auxiliary smooth functions (j)[^ to the true one-electron KS orbitals through a 
linear transformation operator T. As a result, the orthogonality of the original KS orbitals 
is replaced by the orthogonality of the pseudo functions with respect to the overlap operator 
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s = ftr 


), (36) 

and the expectation valne of any local operator is given by 

(6) = (*|6|*) = . (37) 

The KS eigenvalue equation transforms then into a generalized eigenvalue equation 

IcjyD ( 38 ) 

where the transformed Hamiltonian = T^HgT contains extra local and non-local pseu¬ 
dopotential projectors. In the time-dependent case, also the operator needs to be trans¬ 
formed 


= ( 39 ) 

which leads to 

(//r + P) li-P) = l^p) (40) 

The additional operator P = —iT^^ depends on the velocities of the nuclei and for TDDFT 
calculations with immobile ions is identical to zero. When integrating the rt-TDDFT method 
into Ehrenfest dynamics calculations with mobile atoms, explicit calculation of the P oper¬ 
ator is necessary.^l^ 

For the time propagation of Eq. (jlQ]), we use the hrst-order Crank-Nicolson method,^ 
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where the KS orbitals are calculated at each time step from 


*5 + 



+ At) 





(41) 


Since we are using plane wave basis sets, the dimension of Hs{t) and S matrices is usually 
too large to allow a numerically stable inversion and Eq. fl4T|) is solved using the Conjugated 
Gradient Square method. The advantage of the Crank-Nicolson method for the time prop¬ 
agation is the explicit norm conservation of the wave function. Since we are using the hrst 
order version of the method, we forgo the predictor-corrector step and approximate Hs{t+^) 
by hfs(t). Thus, very small time steps must be used and in all calculations presented here 
we have used the 2 attoseconds. 


4 Results and Discussion 

4.1 Optical spectra 

As a hrst application of the novel method, we will consider the calculation of an optical 
spectrum between two coupled chromophores. An optical spectrum can be obtained by 
applying a time dependent electric held to the system such that it samples all frequencies. 
In order to sample all frequencies equivalently, the electric held has ideally the form of a 
5-function, E{t) oc — to). This can be done in two ways: one is adding an electric held in 
the form of a very narrow gaussian. The potential that is added to the periodic Hamiltonian 
of the system and multiplied in real space with a saw function, which drops from 1 to 0 in 
the part of the supercell furthest away from the nuclei. The second way is particularly suited 
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for periodic systems: it prescribes to multiply the occupied KS states by a phas^ 

0i(r, t = 0+) = t = 0") (42) 

where E is the held. For molecular calculations, both methods result in identical time 
evolution of the dipole moment provided that the gaussian pulse integrates to the held 
strength \E\. The dipole moment is calculated at each time step using 

=/p(r,t)rdr (43) 

The oscillator strength is obtained by a Fourier transformation of the time-dependent dipole 
moment into the frequency domain. 

'5(w) = “ Im[akk{,uj)] (44) 

3 TT ^—' 

km 

2 /” _ 2 

Im[akk{uj)] = --pr / sin(a;t)e"'^* [pik{t) - fik{to)]dt (45) 

J 

where k stands for one of the Cartesian directions x,y or z and 7 is a small damping factor.^^Sl 
We choose here a Gaussian damping function rather than the usual exponential as it produces 
cleaner spectra. As a result, the shape of the spectral bands is no longer Lorentzians. In a 
subsystem calculation, the held is applied to each subsystem and the total optical spectrum 
of the system is obtained by summing the oscillator strengths of the subsystems because 

III 

We have calculated the optical spectra of a parallel-stacked benzene fulvene dimer (Fig. 
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d]). The dimers were constructed by placing the relaxed monomers along the Z-axis with 
an intermolecular distance of 4, 5, 6, 7 and 8 A. The PBE functionaP^ and the ultrasoft 
pseudopotentials from the GBRV library^SSl where used throughout. For all calculations, a 
supercell of dimensions 31.0 x 32.5 x 37.8 a.u.^ was used with a kinetic energy cutoff of 55.0 Ry 
and density cutoff of 600.0 Ry. Each rt-TDDFT calculation was run for 16 femtoseconds with 
a time step of 2 attoseconds. The simulation time was tested for convergence by performing 
the Fourier transformation after each femtosecond until no differences in the resulting optical 
spectra were observed. The isolated monomers were calculated using the same supercell and 
cutoff. In the subsystem calculations, the LC94®I functional for the nonadditive kinetic 
energy was used throughout. 

Figure 1: Geometry of the Benzene-Fulvene parallel stacked dimer 


y 
t 



Table [T] lists the interaction energy in the dimer obtained using KS-DFT and subsystem 
DFT for the range of separation distances and the number of electrons misplaced by FDE, 
(Ap) dehned a^HSI 



|Ap(r)| dr 


(47) 


where 

Ap(r) = pFDE(r) - PKs(r) (48) 


KS-DFT predicts a repulsive interaction energy for all distances, consistent with the known 
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Table 1: The interaction energy in the benzene-fulvene dimer calculated using KS-DFT 
and subsystem DFT and the number of misplaced electrons [Eq. flT7|l ] for the different 
intermolecular distances 


R 

(A) 

'^int 

(kcal/mol) 

^int 

(kcal/mol) 

(Ap) 

4 

1.83 

-1.01 

0.016 

5 

0.48 

0.23 

0.002 

6 

0.39 

0.38 

0.000 

7 

0.28 

0.28 

0.000 

8 

0.19 

0.19 

0.000 


dehciency of semilocal GGA functionals to account for van der Waals interactions typical 
for such TT — TT stacked complexes.!^ Subsystem DFT, on the other hand, predicts a negative 
interaction energy at the shortest considered intermolecular distance of 4 A. Such behavior 
of subsystem DFT has been reported and analyzed in detail recent lyBSlI M lI l lllIllMI ^nd while 
the negative interaction energy is more consistent with the physical interactions present in 
the system, it is most likely the result of error compensation of the NAKE which struggles 
with the present density overlap at this short intermolecular distance. For intermolecular 
distances of 6 A and larger, the difference between DFT and subsystem DFT interaction 
energies becomes negligible. This is conhrmed by the (Ap) values which show that for these 
distances, FDE reproduces the KS-DFT density. 

The accuracy of the subsystem rt-TDDFT method, compared to supermolecular rt- 
TDDFT, is depicted in Fig. |2l where the total optical absorption spectrum of the benzene- 
fulvene dimer obtained using the two methods is shown for the three separation distances 
i? = 4, 6 and 8 A. For = 4 A, the difference between the supermolecular and subsystem 
calculation is the largest, which is consistent with the (Ap) value in Table[T]of the dimer from 
the ground state calculations and can be attributed to the failure of the nonadditive kinetic 
energy functionals.^ For the intermolecular distance of i? = 6 A, both methods reproduce 
very similar spectra, with the supermolecular calculation showing three major excitations 
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Figure 2: Optical spectra of the benzene-fulvene dimer at the separation distances of i? = 4, 
6 and 8 A, obtained using rt-TDDFT and subsystem rt-TDDFT methods. 



























































at 4.99 , 6.78 and 7.73 eV and snbsystem calcnlation having very similar valnes of 4.99, 
6.77 and 7.78 eV. These deviations are well within the accnracy of 0.06 eV established for 
the FDE formnlation of the LR-TDDFT method.^ The oscillator strength at the freqnency 
of the lowest/highest excitation energies are identical to KS-DFT and a slight difference is 
fonnd for the middle excitation energy. Snrprisingly, the sitnation becomes worse for the 
largest intermolecnlar distance of 7? = 8 A, where the spectrnm is completely identical be¬ 
low 7.2 eV bnt diverges for the higher freqnencies, where the KS-DFT calcnlation prodnces 
a flattened spectrnm bnt FDE retains the peak at 7.73 eV as in the spectrnm of 77 = 6 
A. In fact, the spectra obtained for the two intermolecnlar distances nsing snbsystem DFT 
are qnasi identical, snggesting that FDE overestimates the intermolecnlar interaction at the 
larger separation distance compared to snpermolecnlar rt-TDDFT. This effect can also be 
seen from the time evolntion of the dipole moment in the snbsystem and snpermolecnlar 
calcnlations, which are depicted for the two intermolecnlar distances in Fig. |3]for the hrst 10 
femtoseconds . Only the x component of the dipole moment is shown, from which the peak 
in qnestion originates. As one can see, there is a discernible difference in the evolntion of 
the dipole moment for the two distances in the snpermolecnlar calcnlation, while in the snb¬ 
system calcnlation the difference between the two distances is very small and both resemble 
the evolntion of the dipole moment in the snpermolecnlar calcnlation for 7? = 6 A. Fnrther 
research is reqnired to determine the physics behind this phenomenon. 

Additional to reprodncing the snpermolecnlar rt-TDDFT resnlt in a divide-and-conqner 
manner, the snbsystem approach has the interesting featnre of allowing to compare properties 
of isolated systems with embedded systems, providing hrst hand information on the effect 
of the embedding on the properties of the snbsystems. Fig. 0] shows the optical spectra of 
the isolated benzene and fnlvene molecnles, compared to their embedded eqnivalents in the 
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dipole strength function (eV'“') dipole strength function (eV’“') 


Figure 3: The time evolution of the dipole moment in the x direction for the intermolecular 
distances of 6 and 8 A, obtained using rt-TDDFT and subsystem rt-TDDFT methods. 
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dimers at different intermolecular separations. The benzene molecule only shows one peak 
at 6.77 eV for the isolated molecule. At = 8 A, the subsystem calculation reproduces 
the supermolecular result exactly, at = 6 A the position of the peak is correct but the 
oscillation strength is slightly underestimated and at 7? = 4 A there is also a slight interaction 
induced shift to the lower frequency of 6.75 eV. The optical spectrum of the fulvene molecule 
is more complicated due to lower symmetry, with four major excitations at 4.99, 6.44, 7.30 
and 7.90 eV. The first three excitations originate from the dipole response to an electric 
held in the x direction, while the last one stems from a response to an electric held in the 
^/-direction. The optical spectra from the subsystem calculation show a very clear ehect 
from the embedding potential generated by the presence of the benzene molecule even at the 
separation distance of 8 A. For the lowest two peaks and the highest peaks, the subsystem 
spectra are shifted to lower energies for the smallest intermolecular distance of i? = 4. The 
shift is smallest for the excitation energy of 4.99 eV (0.05 eV) and largest for the excitation 
energy of 6.44 (0.17). For the intermolecular distances of 7? = 6 and 8, the correct excitation 
energy is reproduced and the ehect of the embedding is visible only in the lower oscillator 
strength, which slowly converges to the value in the isolated molecule. At the highest 
excitation energy of 7.90 eV, there is also an interaction induced shift of 0.13 eV at the 
R = Q and 8 A and the intensity is overestimated for the shorter distances. The situation 
for the peak at 7.30 eV is diherent, where all embedded fulvene spectra show similar shift of 
0.13 eV (0.14 eV for R = 4). Additionally, the oscillator strength diverges from the isolated 
value with increasing intermolecular distance. 

The question of the interaction induced shift has been previous discussed within FDE 
using the linear response formulation of the method jSSl I^rni lSlin] the original implementa¬ 
tion of linear response FDE,^ the linear density response of a subsystem is obtained while 
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Figure 4: Optical spectra of the isolated benzene and fulvene molecules, obtained using rt- 
TDDFT, compared to the spectra of the embedded molecules in the dimers at the separation 
distances of i? = 4, 6 and 8 A., obtained using subsystem rt-TDDFT. 
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keeping the densities of the other subsystems frozen in time. As a result, the response is re¬ 
stricted to the active subsystem only.®®] As mentioned above, this uncoupled response, Xj, 
is equivalent to performing a subsystem rt-TDDFT calculation where only the KS orbitals of 
one subsystems are propagated in time, while the other subsystems are kept in their initial 
state. The advantage of such an approach is the considerable computational simplicity. Nu¬ 
merous applications on optical spectra,®^ Raman spectra!^ induced circular dichroisnJ^I^ 
and electron-spin-resonance hyperhne coupling^ have shown that the uncoupled response is 
sufficient in reproducing supermolecular results with good accuracy, even in the presence of 
hydrogen bonds, as long as there are no couplings in the excitation between the systems. An 
example of when the uncoupled response is not enough is the case of identical chromophores, 
where explicit excitation transfer occurs, as will be discussed in Section 021 In later work, J. 
Neugebauei® introduced a framework for including coupling between selected excitations in 
the linear response regime. The subsystem rt-TDDFT method offers a straightforward way 
to include couplings between all excitations, as well as reproducing the uncoupled results. 
Clearly, a fully coupled calculation is more computationally expensive than an uncoupled 
one. 

The difference in the results between the coupled and uncoupled calculations for the 
benzene-fulvene dimer is shown in Figs. [5] and [6] for different intermolecular separations. 
The largest differences between the coupled and uncoupled results are found for the fulvene 
molecule for the excitation energy at 6.44 eV for both the interaction induced shift and the 
oscillation strength. In the fully coupled calculation, the dynamical embedding potential due 
to the presence of the benzene molecule results in a shift of 0.20 eV at the separation distance 
of 4 A, while no shift is present for the uncoupled calculation (i.e. stationary embedding). 
Additionally, a large difference in the oscillation strength for this excitation is found between 
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Figure 5: Optical spectra of the embedded benzene molecule, obtained using coupled and 
uncoupled subsystem rt-TDDFT methods at the separation distances of i? = 8, 6 and 4 A. 
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Figure 6: Optical spectra of the embedded fulvene molecule, obtained using coupled and 
uncoupled subsystem rt-TDDFT methods at the separation distances of i? = 8, 6 and 4 A. 
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the coupled and uncoupled calculations. Though not as pronounced, a similar effect is found 
for the 4.99 eV excitation in fulvene and 6.78 eV in benzene. At = 4 A, the coupled 
calculation produces an interaction induced shift of 0.04 eV, compared to 0.02 eV in the 
uncoupled calculation, and a lower oscillation strength. At the intermediate distance of 6 
A, the effect of the lower oscillation strength is still present, but no interaction induced 
shift is observed. At 8 A, the coupled and uncoupled spectra are almost identical, with the 
exception of the 6.44 eV excitation in the fulvene molecule, where there is still a small effect 
on the oscillation strength. It should also be noted that for the uncoupled excitations in 
the fulvene spectrum, the coupled calculation produces slightly larger oscillation strengths, 
which is necessary in order to satisfy the sum rule of the spectrum. 

4.2 Excitation energy transfer 

In the last section, we have seen that subsystem rt-TDDFT allows us to reproduce full cou¬ 
pling between subsystems in a straightforward fashion. A strongly coupled system which is 
often studied using TDDFT is a dimer consisting of two identical chromophores. An ade¬ 
quate description of such system requires a method which includes the full coupled response 
between all subsystem,^ as clearly shown by Eq. fl34|) .^ The coupling in excitations along a 
certain direction between the chromophores opens a window for an explicit excitation energy 
transfer (EET). The central quantity in EET is the electronic coupling V, defined as 

V = (D*AIHIDA*) (49) 

where \D*A) and \DA*) are orthogonal diabatic states. Herein, D and A represent the donor 
and acceptor in the excitation energy exchange transfer process and H is the molecular 
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Hamiltonian operator.l^l^^ The coupling is linked to the excitation energy transfer rate 
through Fermi’s Golden Ruld^ 

= 2ti\V\^ j J{e)de (50) 

where J(e) is the spectral overlap between the absorption spectra of the D and A chro- 
mophores. Within the framework of linear response TDDFT, the electronic coupling can be 
obtained from the Davydov splittin^^l^ of the coupled excitation energy in the spectrum, 
compared to the absorption spectrum of the monomer. Real time TDDFT, however, offers 
the advantage to study such transport events in real time and this can offer an intuitive pic¬ 
ture of the EFT process. For this purpose, we chose the test system of two sodium dimers, 
situated along the z-axis with separation distances of 12, 15 and 17.5 Bohr. This system 
has been well studied beford^^^l^ and is known to couple strongly at the excitation energy 
of 2.18 eV in the direction of the Na-Na bond. The full cluster was placed in a supercell 
of 22.7 X 22.7 x 43.5 a.u.^ to avoid periodic interactions. The calculation was performed 
with ultrasoft pseudopotentials from the GBRV® library with a kinetic energy cutoff of 55.0 
Ry and density cutoff of 660.0 Ry and the LGO^®^ functional was used for the nonadditive 
kinetic energy. The exact value of the excitation energy was hrst obtained from an opti¬ 
cal absorption spectrum of the Na 4 cluster using FDE. At the start of the simulation, the 
\D*A) state is prepared by applying an electric held in a direction along the Na-Na bond to 
only one of the subsystems, noted as the donor, with a frequency of 2.18 eV, using a cosine 
function damped by a Gaussian. This state is not stationary and will evolve accordingly by 
populating the \DA*) state. The time evolution of the held is depicted in Fig. [3 Both of 
the subsystems are evolved simultaneously for 100 fs with a time step of 2 as. At each time 
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step, the dipole moment in the direction of the excitation is recorded for both the donor and 
acceptor molecule, which are depicted in Fig. [8] for the three separation distances. 

Figure 7: Applied electric held to the Na 4 cluster. 



Starting from f = 0, the dipole moment of the donor molecule reacts to the applied held 
and the dipole moment starts oscillating, reaching a maximum within the 2.5 fs, at which 
point the \D*A) state reaches maximum population. From that point on, the system proceeds 
in transferring the excitation energy acquired from the applied held to the \DA*) state. It 
is important to note that the acceptor subsystem has no direct held applied in its potential 
and the excitation transfer occurs solely through the response of the embedding potential to 
the perturbation. At the shortest separation distance of 12 Bohr, the excitation energy is 
fully transfered to the \DA*) state within approximately 20 fs, 13 fs after the damping of the 
applied held has taken place. The complete transfer of the excitation energy is clear from 
the intensity of the dipole moment when it reaches its maximum values, which is equal to the 
maximum in the population of \D*A) and \DA*) states. As no nuclear relaxation is included 
in this simulation, the excitation energy is transfered back to the \D*A) state. As can be seen 
from Fig. [8], the excitation rate strongly depends on the separation distance. According to 
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Figure 8: Excitation energy transfer between two Na 2 molecules at the separation distances 
of 12.5, 15 and 17.5 Bohr, after the excitation of the donor molecule with an applied electric 
held of 2.18 eV frequency. 
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Forster theoryj^^the excitation transfer rate has a dependence with n = 6 at long range 
and a lower n dependence for shorter distances due to exchange interactions, as described by 
Dexter theory.^ This has been previously conhrmed for this particular system by Hofman 
et ah,^ which found the Forster theory to be valid for separation distances of 25 Bohr and 
higher. This coincides with the value of n = 5.1 observed here for the shorter separation 
distances. The excitation energy transfer rates extracted from Fig. |8] are summarized in 
Table [2J 

Table 2: Excitation energy transfer (EET) rate for the Na 4 cluster at the different separation 
distances. 


R 


(Bohr) 

(fs)-i 

12.5 

0.052 

15.0 

0.022 

17.5 

0.009 


5 Summary and Conclusions 

We presented an extension of the subsystem DFT method to real-time TDDFT. In the new 
method, the electron density is evolved in time alongside a time dependent potential as a 
collection of time-dependent subsystem electron densities. The implementation is such that 
one can choose between evolving only one of the subsystems in time while keeping the other 
subsystems frozen at the initial time to (uncoupled simulation) or evolving all subsystems 
simultaneously (coupled simulation). The former choice will only include the response of the 
excitations within the subsystem, while the latter will include all the couplings between all 
of the subsystems. 

The method was first applied to obtain optical absorption spectra of a benzene-fulvene 
dimer, where similar accuracy w.r.t. supramolecular DFT calculations were found as for 
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the case of ground state FDE and LR-TDDFT FDE applications. For shorter distances, 
the optical absorption spectra start to diverge from the supramolecular calculations due 
to the known failure of NAKE functionals, while for intermediate and long distances the 
supramolecular results are reproduced with good accuracy. 

The choice between performing a coupled and an uncoupled calculation presented the 
opportunity to examine the effect of the excitation coupling between the subsystems on 
the reproduction of optical absorption spectra. A clear and consistent effect is seen on 
the interaction induced shift and the oscillator strength when omitting the full coupling of 
the subsystems, for the coupled excitations. Since we only examined a pair of nonidentical 
chromophores here, this effect is limited and we were able to conhrm the previous Ending^ 
that subsystem LR-TDDFT is able to reproduce the supramolecular results for most systems 
while neglecting the intersubsystem coupling. 

As a second application, the real-time excitation energy transfer was studied in a Na 4 
cluster, where only one of the Na 2 was subjected to a periodic electric held with a frequency 
of a coupled excitation energy. We found that the embedding potential is fully capable of 
transmitting the excitation energy between the subsystems and that the rate of the excitation 
energy transfer coincided with the Forster theory prediction. 

The ability to represent the real-time dynamics of excitation energy transfer (EET) be¬ 
tween subsystems is an important Ending of this work. The ability to generate diabatic 
states for EET in a way that resembles most the true experiments (i.e., via an applied held) 
is of pivotal importance, complementing the linear-response theories. Future directions to¬ 
wards extending the method to period systems and to the inclusion of nuclear dynamics 
(Ehrenfest) are undergoing. 
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